
*
*adaptation of De Loecker's code
*De Loecker, J. 2013. Detecting Learning by Exporting, American Economic Journal: Microeconomics.
*-Code last saved on December 10, 2012
*-----------------------------------------------------------------------------------------------------*

***
*** COMMENT: SOMETIMES THE CODE DOES NOT RUN. IT PRODUCES 
*** ERROR r(3000) "variable undeclared" or error r(3499) "MODEL1() not found"
*** However, it is a matter of closing stata and opening it again for it to
*** run properly
*** 
***

*upload
clear
clear matrix
clear mata
set more off
use $analysis_data , clear 
drop y 
rename va y

*variables' names and creation
gen Id = ID
destring Id, replace
gen exportdum = (year>=event)
rename sector_g nace

*Creating event dummies
sort Id year
gen years_export = year if exportdum==1

*******************************************
*it is better to use the lag fuction than just _n-1 since it mix values
*for different firms 
*******************************************
xtset Id year

*we assume a 0.2 depreciation rate
gen i = k - 0.8*l.k
gen exportdum_lag = l.exportdum
keep year Id l k i m y exportdum nace D_* event

*******************************************
*keeping non empty variables
keep if y!=.
*******************************************

*note: the i.x*i.z option was not working
*note: the article states that the first reg of k and i was supposed to include l (Table 1 note)
gen year_nace = year*nace
reg k exportdum l i.year_nace
reg k exportdum l y i.year_nace
reg i exportdum l i.year_nace
reg i exportdum l y i.year_nace
*Panel B
* set the panel accordingly:
xtset Id year, yearly
forvalues i=1/4 {
gen dk`i'=k-l`i'.k
reg dk`i' exportdum
}
*-----------------------------------------------------------------------------------------------------------*
* Section III: Results
*1. Estimate Production functions by industry (nace)
*2. Estimate DID 
*3. Estimate Nonparametric
*4. Estimate Extension g(e,omega,i) 
*---------------creating variables-------------------------------------------------------------------*
* higher order terms on inputs
local M=3
local N=3
forvalues i=1/`M' {
gen l`i'=l^(`i')
gen m`i'=m^(`i')
gen k`i'=k^(`i')
*interaction terms
forvalues j=1/`N' {
gen l`i'm`j'=l^(`i')*m^(`j')
gen l`i'k`j'=l^(`i')*k^(`j')
gen k`i'm`j'=k^(`i')*m^(`j')
}
}
gen lkm=l*k*m
gen l2k2m2=l2*k2*m2
gen l3k3m3=l3*k3*m3
*----------------------------------------------------------------------------------------------------*
* FIRST STAGE ACF- with materials 
*[note 3rd order polynomial approximation is used to keep comparable to e.g. Levinsohn and Petrin (2003) and levpet.do]
xi: reg y l* m* k* exportdum i.year
predict phi
gen phi_lag=L.phi
*----------------------------------------------------------------------------------------------------*
gen l_lag=L.l
gen k_lag=L.k
gen l_lag2=l_lag^2
gen k_lag2=k_lag^2
gen l_lagk_lag=l_lag*k_lag
gen lk=l*k
gen l_lagk=l_lag*k
gen inv_exp=i*exportdum
gen exp_lag= L.exportdum
gen i_lag= L.i


drop _I*
sort Id year
gen constant=1
drop if y==.
drop if l==.
drop if k==.
drop if i==.
drop if phi==.
drop if k_lag==.
drop if i_lag==.
drop if l_lag==.
drop if phi_lag==.

/* Overview specifications:
1. static control (m) and exogenous process
2. static control (m) and endogenous process (DL1): DL1
3. static control (m) and endogenous process (DL2) with action function (Eq 10) section III: DL2

* Note flexibility regarding inclusion of (lagged) export participation, here I just show the use of lagged export dummy.
Results in paper are robust to the inclusion of additional variables such as e.g. export intensity.
*/
*-------------------------------------BEGIN MATA PROGRAM--------------------------------------------*
mata
mata set matastrict off
void GMM_1(todo,betas,crit,g,H)
{
	PHI=st_data(.,("phi"))
	PHI_LAG=st_data(.,("phi_lag"))
	Z=st_data(.,("constant","l_lag","k"))
	X=st_data(.,("constant","l","k"))
	X_lag=st_data(.,("constant","l_lag","k_lag"))
	Y=st_data(.,("y"))
	C=st_data(.,("constant"))
	
	OMEGA=PHI-X*betas'
	OMEGA_lag=PHI_LAG-X_lag*betas'
    OMEGA_lag2=OMEGA_lag:*OMEGA_lag
    OMEGA_lag3=OMEGA_lag2:*OMEGA_lag
    OMEGA_lag_pol=(C,OMEGA_lag,OMEGA_lag2,OMEGA_lag3)
	g_b = invsym(OMEGA_lag_pol'OMEGA_lag_pol)*OMEGA_lag_pol'OMEGA
	XI=OMEGA-OMEGA_lag_pol*g_b
	crit=(Z'XI)'(Z'XI)
}

void GMM_DL1(todo,betas,crit,g,H)
{
	PHI=st_data(.,("phi"))
	PHI_LAG=st_data(.,("phi_lag"))
	Z=st_data(.,("constant","l_lag","k"))
	X=st_data(.,("constant","l","k"))
	X_lag=st_data(.,("constant","l_lag","k_lag"))
	Y=st_data(.,("y"))
	C=st_data(.,("constant"))
	EXP=st_data(.,("exportdum"))
	EXP_lag=st_data(.,("exp_lag"))

    
	OMEGA=PHI-X*betas'
	OMEGA_lag=PHI_LAG-X_lag*betas'
    OMEGA_lag_EXP=OMEGA_lag:*EXP_lag
    OMEGA_lag2=OMEGA_lag:*OMEGA_lag
    OMEGA_lag3=OMEGA_lag2:*OMEGA_lag
    OMEGA_lag2_EXP=OMEGA_lag2:*EXP_lag
    OMEGA_lag3_EXP=OMEGA_lag3:*EXP_lag
    OMEGA_lag_pol=(C,OMEGA_lag,OMEGA_lag2,OMEGA_lag3,OMEGA_lag_EXP, OMEGA_lag2_EXP,OMEGA_lag3_EXP,EXP_lag)
	g_b = invsym(OMEGA_lag_pol'OMEGA_lag_pol)*OMEGA_lag_pol'OMEGA
	XI=OMEGA-OMEGA_lag_pol*g_b
	crit=(Z'XI)'(Z'XI)
}

void GMM_DL2(todo,betas,crit,g,H)
{
	PHI=st_data(.,("phi"))
	PHI_LAG=st_data(.,("phi_lag"))
	Z=st_data(.,("constant","l_lag","k"))
	X=st_data(.,("constant","l","k"))
	X_lag=st_data(.,("constant","l_lag","k_lag"))
	Y=st_data(.,("y"))
	C=st_data(.,("constant"))
	EXP=st_data(.,("exportdum"))
	EXP_lag=st_data(.,("exp_lag"))
	INV_lag=st_data(.,("i_lag"))

	OMEGA=PHI-X*betas'
	OMEGA_lag=PHI_LAG-X_lag*betas'
	OMEGA_lag2=OMEGA_lag:*OMEGA_lag
    OMEGA_lag3=OMEGA_lag2:*OMEGA_lag
    OMEGA_lag_INV=OMEGA_lag:*INV_lag
    OMEGA_lag_EXP=OMEGA_lag:*EXP_lag
    OMEGA_lag_IE=OMEGA_lag_INV:*EXP_lag
	INV_EXP_lag=INV_lag:*EXP_lag
	OMEGA_lag_pol=(C,OMEGA_lag,OMEGA_lag2,OMEGA_lag3,OMEGA_lag_INV,OMEGA_lag_EXP,OMEGA_lag_IE,INV_lag,EXP_lag,INV_EXP_lag)
	g_b = invsym(OMEGA_lag_pol'OMEGA_lag_pol)*OMEGA_lag_pol'OMEGA
	XI=OMEGA-OMEGA_lag_pol*g_b
	crit=(Z'XI)'(Z'XI)
}

void GMM_DL_linear(todo,betas,crit,g,H)
{
	PHI=st_data(.,("phi"))
	PHI_LAG=st_data(.,("phi_lag"))
	Z=st_data(.,("constant","l_lag","k"))
	X=st_data(.,("constant","l","k"))
	X_lag=st_data(.,("const","l_lag","k_lag"))
	Y=st_data(.,("y"))
	C=st_data(.,("constant"))
	EXP=st_data(.,("exportdum"))
	EXP_lag=st_data(.,("exp_lag"))

	OMEGA=PHI-X*betas'
	OMEGA_lag=PHI_LAG-X_lag*betas'
	OMEGA_lag_pol=(C,OMEGA_lag,EXP_lag)
	g_b = invsym(OMEGA_lag_pol'OMEGA_lag_pol)*OMEGA_lag_pol'OMEGA
	XI=OMEGA-OMEGA_lag_pol*g_b
	crit=(Z'XI)'(Z'XI)
}


void MODEL1()
	{
S=optimize_init()
optimize_init_evaluator(S, &GMM_1())
optimize_init_evaluatortype(S,"d0")
optimize_init_technique(S, "nm")
optimize_init_nmsimplexdeltas(S, 0.1)
optimize_init_which(S,"min")
optimize_init_params(S,(2,0.8,0.2))
p=optimize(S)
p
st_matrix("beta1",p)
}


void MODEL_DL1()
	{
S=optimize_init()
optimize_init_evaluator(S, &GMM_DL1())
optimize_init_evaluatortype(S,"d0")
optimize_init_technique(S, "nm")
optimize_init_nmsimplexdeltas(S, 0.1)
optimize_init_which(S,"min")
optimize_init_params(S,(2,0.8,0.2))
p=optimize(S)
p
st_matrix("beta_NP",p)
}

void MODEL_DL2()
	{
S=optimize_init()
optimize_init_evaluator(S, &GMM_DL2())
optimize_init_evaluatortype(S,"d0")
optimize_init_technique(S, "nm")
optimize_init_nmsimplexdeltas(S, 0.1)
optimize_init_which(S,"min")
optimize_init_params(S,(2,0.8,0.2))
p=optimize(S)
p
st_matrix("beta3",p)
}

void MODEL_DL_linear()
	{
S=optimize_init()
optimize_init_evaluator(S, &GMM_DL_linear())
optimize_init_evaluatortype(S,"d0")
optimize_init_technique(S, "nm")
optimize_init_nmsimplexdeltas(S, 0.1)
optimize_init_which(S,"min")
optimize_init_params(S,(2,0.8,0.2))
p=optimize(S)
p
st_matrix("beta_lin",p)
}

end
* note that in the past optimizations the starting values are simply entered. Use OLS estimates as starting values.

//produce parameters on law of motion, directly related to decomposition of total productivity effect
*---------------------------------------------------------------------------------------------------------------*
sort Id year
mata MODEL1()
gen const=beta1[1,1]
gen betal=beta1[1,2]
gen betak=beta1[1,3]
gen omega=phi-betal*l-betak*k

mata MODEL_DL1()
gen const_NP=beta_NP[1,1]
gen betal_NP=beta_NP[1,2]
gen betak_NP=beta_NP[1,3]
gen omega_NP=phi-betal_NP*l-betak_NP*k


gen omega2_NP=omega_NP^2
gen omega3_NP=omega_NP^3
gen omega_NP_exp=omega_NP*exportdum
gen omega2_NP_exp=omega2_NP*exportdum
gen omega3_NP_exp=omega3_NP*exportdum

reg omega_NP L.omega_NP L.omega2_NP L.omega3_NP L.exportdum L.omega_NP_exp L.omega2_NP_exp L.omega3_NP_exp
gen LBE_NP=_b[L.exportdum]+_b[L.omega_NP_exp]*omega_NP+_b[L.omega2_NP_exp]*omega2_NP+_b[L.omega3_NP_exp]*omega3_NP

*----------------------------------------------------------------------------------------------------------------*
mata MODEL_DL2()
gen const3=beta3[1,1]
gen betal3=beta3[1,2]
gen betak3=beta3[1,3]
gen omega3=phi-betal3*l-betak3*k-const3
gen omega3_2=omega3^2
gen omega3_3=omega3^3
gen int1=exportdum*i
gen int2=omega3*i
gen int3=omega3*exportdum
gen int4=omega3*i*exportdum
reg omega3 L.omega3 L.omega3_2 L.omega3_3 L.i L.exportdum L.inv_exp L.int1 L.int2 L.int3 L.int4
gen theta3_1=_b[L.omega3]
gen theta3_2=_b[L.omega3_2]
gen theta3_3=_b[L.omega3_3]
gen theta3_4=_b[L.i]
gen theta3_5=_b[L.exportdum]
gen theta3_6=_b[L.int1]
gen theta3_7=_b[L.int2]
gen theta3_8=_b[L.int3]
gen theta3_9=_b[L.int4]

*merging back 
tostring Id, gen(ID) format(%24.0g)
gen prod_DL = omega_NP
keep ID year prod_DL 
tempfile prod
save `prod', replace

use $analysis_data , clear 
merge 1:1 year ID using `prod'
drop _m

save temp/temp_DL.dta , replace
